Quantum Monte Carlo Calculations of Neutron Matter 



J. Carlson 

Theoretical Division, 
Los Alamos National Laboratory, 
Los Alamos, New Mexico 87545 

J. Morales, Jr, V. R. Pandharipande and D. G. Ravenhall 

Department of Physics, University of Illinois at Urbana- Champaign, 
1110 W. Green St., Urbana, IL 61801, U.S.A. 
(Dated: February 8, 2008) 

Abstract 

Uniform neutron matter is approximated by a cubical box containing a finite number of neutrons, 
with periodic boundary conditions. We report variational and Green's function Monte Carlo calcu- 
lations of the ground state of fourteen neutrons in a periodic box using the Argonne t;8' two-nucleon 
interaction at densities up to one and half times the nuclear matter density. The effects of the finite 
box size are estimated using variational wave functions together with cluster expansion and chain 
summation techniques. They are small at subnuclear densities. We discuss the expansion of the 
energy of low-density neutron gas in powers of its Fermi momentum. This expansion is strongly 
modified by the large nn scattering length, and does not begin with the Fermi-gas kinetic energy 
as assumed in both Skyrme and relativistic mean field theories. The leading term of neutron gas 
energy is ~ half the Fermi-gas kinetic energy. The quantum Monte Carlo results are also used to 
calibrate the accuracy of variational calculations employing Fermi hypernetted and single operator 
chain summation methods to study nucleon matter over a larger density range, with more realistic 
Hamiltonians including three-nucleon interactions. 

PACS: 21.65.+f, 26.60.+C 
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I. INTRODUCTION 



Since the discovery of neutron stars in 1967 there has been a continued interest in cal- 
culating the properties of neutron matter from realistic models of nuclear forces 0, Q]- It 
is difficult to extrapolate the data on bound nuclei using energy-density functionals to es- 
timate the equation of state E{p) of pure neutron matter. Various Skyrme and relativistic 
energy-density functionals which fit the binding energies and radii of nuclei available in 
laboratories at present give rather different E{p) for neutron matter. These energy-density 
functionals also predict different properties of nuclei near the neutron drip line |2,], which 
may be synthesized in the near future using radioactive ion beams. Theoretical predictions 
of the neutron matter E{p) have been used to constrain the energy density functionals used 
to study neutron rich nuclei. 

The two-neutron interaction is strong and highly spin dependent. Therefore calculating 
the neutron matter E{p) is a challenging many-body problem, though in some ways it 
is simpler than that of symmetric nuclear matter. Neutron matter properties have been 
calculated recently with Brueckner theory 0, |^ and with variational methods using chain 
summation techniques [1, 0]. There is good agreement between the results of these two 
methods and recent high precision models of nn interaction give rather similar neutron 
matter E(p) with the lowest order Brueckner method 0|. The results for symmetric nuclear 
matter, however, have more model dependence. The E{p) of high density neutron matter is 
also sensitive to the lesser known three neutron interaction 0. 

The Brueckner and variational methods use different expansions. From the results of 
two- and three-hole line contributions, the cluster expansion in Brueckner theory appears to 
converge though contributions of clusters with more than three neutrons have not been 
calculated. In contrast, the cluster expansion of the energy expectation value of neutron 
matter used in the variational method has a rather poor convergence when the optimum 
variational wave function, is used. In the latest calculations two- and three-body 
cluster contributions are calculated accurately, while those of > 4-body clusters are summed 
approximately with hypernetted and single-operator chain summation methods. The con- 
vergence rate of the expansion is sensitive to the range of correlations in "^y. Hence it is often 
possible to use shorter range correlations, which give a more convergent cluster expansion 
together with a variational energy within a few % of the optimum minimum. Thus, even 
though the results of these two methods are in agreement within several %, the theoretical 
error in the treatment of long range correlations is not well estimated. 

In the past few years it has been possible to calculate the energies of all the bound states 
of nuclei having up to ten nucleons with errors estimated to be ^ 2 % using the Green's 
function Monte Carlo (GFMC) method IH llll . Results of these calculations are being used 
to construct realistic models of three-nucleon interactions fl3\. The computational effort 
necessary for a nuclear GFMC calculation scales approximately with 2^A\/{N\Z\) for a 
system with N neutrons, Z protons and A = N + Z . The factor 2^ comes from the number 
of spin states of A nucleons and A\ / {N\Z\) is the number of charge conserving isospin states. 
In the present work we report calculations of the ground state of 14 neutrons in a periodic 
box (PB) with the GFMC method considering all the 2^^ spin states. We have also used the 
auxiliary field diffusion Monte Carlo (AFDMC) method proposed by Schmidt and Fantoni 
jisj l to calculate the ground state energy. In this method one effectively samples the 2^^ spin 
states stochastically. The computational effort of AFDMC scales with A^, and thus it can 
be used to study systems with larger values of A 3|. The present AFDMC calculations 
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seem to have larger errors than the GFMC; however, it may be possible to improve their 
accuracy. 

The interactions and the variational wave functions used in this work are described in 
sect. II. The quantum Monte Carlo calculations, variational (VMC), GFMC and AFDMC, 
are described briefly in section III, where we present results at p = 0.04, 0.08, 0.16 and 0.24 
fm~^. The details of these methods have been presented previously; here we simply describe 
the additional techniques used to calculate results for 14 neutrons in a PB, and discuss 
several tests of the calculations. The total energy and the potential energy expectation 
values are reported for each density. 

The variational calculations using chain summation methods (VCS) are reported in sec- 
tion IV. In this section we also discuss the difference between the density matrices of 14 
neutrons in a PB and of uniform gas (UG) with large number of neutrons. The smallness 
of this difference makes the 14-particle PB a useful approximation to UG. The difference 
between the energy per neutron in the PB and in UG is estimated using variational wave 
functions. It is small at subnuclear densities, p < Po = 0.16 fm~^, but significant at 1.5 po. 
The comparison of VCS results with the QMC suggests that the former can have errors up 
to ~ 10 %. 

The pair distribution functions obtained from VMC and GFMC calculations are compared 
in section V. These indicate that neutron matter has strong correlations even at small 
densities, as expected from the large scattering length, a ~ —18 fm in the ^Sq state. The 
results for the E{p), extrapolated to the UG limit, are presented in section VI. Here we 
also discuss the expansion of the E{p) in powers of kp. When lakpl < 1, this expansion 
begins with the Fermi-gas kinetic energy, TpQ = 0.3A;^/m. However, at densities of interest 
in nuclear or neutron star physics {akpl » 1, and the expansion of the E{p) seems to 
begin with ~ Tpc/'^i which is the estimated UG energy for a short range interaction with 
scattering length (—a) oo. This approximation to nuclear forces in low-density neutron- 
gas was suggested by Bertsch [Ts^ . 

The accuracy of the present calculation is discussed in the last section VII. At densities 
< Po the GFMC calculation appears to be well converged and presumably has an accuracy 
of ~ 2 % for the energy of normal neutron matter. However, there is exceptionally strong 
pairing in dilute Fermi gases with akp — oo and their superfluid state can have 
energies below those of the normal state by ~ 10 The main error in the predicted 

energy of low density neutron matter, where three-body forces are small, is likely to be due 
to the neglect of the superfluidity in the present calculation. 



II. INTERACTIONS AND VARIATIONAL WAVE FUNCTIONS 

We have used the Argonne v8' two-nucleon interaction in this work. This simplified 
interaction equals the isoscalar part of the realistic Argonne vl8 interaction in all the S- and 
P- waves as well as in ^Di and its coupling to ^5*1. In neutron matter this interaction can be 
written as an operator with four terms: 

Vij = Mrrj)Ofj , (2.1) 

p=l,4 

Op'^'*-' = 1, ■ tr,, L ■ S . (2.2) 

Here Sij and L • S are the usual tensor and spin-orbit operators. In the calculations using 
the PB boundary condition, the interaction is truncated at Vij = L/2, where L is the length 
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of the cubic box holding 14 neutrons {Li^p = 14) 



v{rij) = v{rij) - Uj) + v{rij) 9{rij - ^) . (2.3) 

The contribution of the long range part, f (rjj) 9{rij — -1), to the E{p) of UG is estimated 
using variational calculations. While at low densities this contribution is small, at p = 0.24 
fm~^ it becomes comparable to the total E{p) in magnitude. 

The variational wave function \&v' used in this work has the form: 



<ilv= \ Sl[F,A<!>, (2.4) 




where $ is the noninteracting Fermion wave function. In UG calculations $ = $_fG) the 
Fermi-gas wave function, while in calculations using the PB boundary conditions $ = $p_b. 
It has 14 neutrons occupying spin up and down states with momenta 

k = 0, ± /cfiX, ± ksy, ± fc^z . (2.5) 

Here ks = 2tt/L and x, y and z are unit vectors. 

The S n denotes a symmetrized product of the noncommuting Fij pair correlation oper- 
ators. In VCS calculations they have four terms involving the four operators of Eq. ()2.2|1 : 

F^,= E Pp fpir.,)Of^ . (2.6) 

p=c,a-,t,b 

The correlation functions fpivij) are obtained by solving two-body Schrodinger-like equations 
and have three parameters, d^dt and a. They correspond to the range of all but the 
tensor correlations (rf), range of tensor correlations {dt), and the average quenching of spin- 
dependent interactions in matter (a). 

In the case of UG, the values of d, dt and a are determined by minimizing the energy with 
the VCS method for (5p = 1. Constraints imposing conservation of nucleons are used during 
this minimization to prevent the Fij from going into regions where the chain summation 
approximation fails. The parameters (3a ■, A and I3h provide additional variation of F; is 
not a variational parameter since F{rij — > oo) = 1. The (3p^c parameters were not used in 
recent calculations [3| since they do not lower the energy significantly after optimizing d, dt 
and a. They are used here for the following reason. The optimum values of d and df in UG 
are > L/2. However, in VMC as well as VCS calculations using the PB boundary condition 
the d and dt must be < L/2. In all PB calculations we use d = dt = L/2, and vary the a 
and jSp^c to minimize the energy. 

In VMC calculations the spin-orbit correlations in the Fij are neglected due to computa- 
tional difficulties associated with the gradient operator in L. These calculations use the v6 
interaction obtained by dropping the spin-orbit term in the v8'. There results are compared 
with those of VCS with the same F to test the accuracy of the chain summation approxima- 
tion. The complete v8' interaction is used in the GFMC calculations where the spin-orbit 
correlations are generated by propagation in imaginary time as discussed in the next section. 
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III. QUANTUM MONTE CARLO CALCULATIONS 



Quantum Monte Carlo methods have often been used to study infinite systems of either 
fermions or bosons at both zero and finite-temperature. Examples include atomic liquid 
^He and ^He [l^, the electron gas, [i^ as well as a myriad studies of lattice models in 
condensed matter theory. They have proven remarkably successful at studying the equation 
of state of strongly-interacting systems, and have also been used to explore phase transitions, 
momentum distributions, static and dynamic response, etc. Although studies of fermion 
systems are usually treated via approximate fixed-node j2l| or constrained path (CP) [i^ 
methods, these approximations can often be quite accurate. 

The nuclear many-body problem is more difficult than all the cases listed above, because 
of the strong spin-isospin dependence of the interaction. Instead of a single function of 
the 3A coordinates of the particles, the wave function of simple systems, the nuclear state 
is described with a set of (complex) amplitudes dependent upon the spins and isospins of 
the nucleons. This complexity has been handled successfully for few-body (A < 10) nuclei 



ll[ by simply summing explicitly over all these amplitudes. Monte Carlo is then used to 
evaluate the 3A-dimensional spatial integrals. 

Variational Monte Carlo (VMC) calculations evaluate the energy and other observables 
thr oug h the use of the Metropolis Monte Carlo method. The method is described in detail 
in 10], the basic idea being to generate points in the 3A-dimensional configuration space 
distributed with the probability density of a weight function iy(R). Here R is the 3A- 
dimensional configuration vector ri, ...r^. The expectation values of operators are obtained 
as averages over the sampled points Rj: 

_ S.(^(R.)|Q|^(R.))/iy(R.) 

The optimum weight function in most cases is the square of the wave function (\E'(R)|\E'(R)) 
for which the denominator of (O) has zero variance. For a system of 14 neutrons the 
(\1'(R)|\E'(R)) is a sum of squares of the 2^'^ spin amplitudes. 

This method grows exponentially in computational time with increasing A, and present- 
day computers limit practical simulations to roughly 14 neutrons. This is somewhat larger 
than the largest nuclei handled to date because there is only one isospin component to the 
wave function. Another limitation of these initial calculations is that we have dropped the 
L ■ S pair correlation functions, as they depend upon the momentum of the particles in 
the pair. A complete evaluation of these terms would be difficult because the derivative 
operators in one pair correlation function can, in principle, act on all other pair correlations. 
This limitation is not very important at low densities, but can be quite significant at higher 
densities. The variational wave function cannot adequately describe p-wave pairing of the 
neutrons which appears to be important at nuclear densities and above. It may be possible 
to construct a simplified wave function which includes most of these correlations in the 
future. 

Green's function Monte Carlo (GFMC) methods are then used to obtain the ground- 
state energy and other properties for the 14 neutrons with periodic boundary conditions. 
The method is the same as that used for light nuclei l3|, with only very minor modifications 



used to implement the periodic boundary conditions. The basic idea is to sample a wave 
function \I'(t) by evaluating path integrals of the form: 

^(r) =llexp[-{H - £;o)Ar]|^y >, (3.2) 
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where each step in the product evolves the system over a short imaginary time At; after 
many steps n — cx), \l/(r = nAr) will converge to the true ground state of the system 
as long as the original (variational) wave function is not orthogonal to it. Because we are 
studying systems at densities higher than equilibrium nuclear density po, the time step used 
here is 0.00025 MeV~\ or 1/2 the time step typically used in nuclear calculations. Again 
the spatial integrals are done with Monte Carlo, using a sum over many configurations with 
different spatial coordinates R. Each configuration includes amplitudes for all the 2^'^ spin 
states which are explicitly summed in the evaluation of matrix elements. 

Some of the GFMC results reported here are obtained with the CP approximation |lol . 2^ . 



Since the neutrons are fermions, they may exchange and produce contributions to the wave 
function of opposite signs, and indeed with arbitrary complex phases. In an exact GFMC 
calculation, this leads to a statistical error that grows with r. This problem is more severe 
at higher densities (or with larger numbers of particles) since it is then easier for a pair to 
interchange. 

To deal with this problem, we implement a constraint on the paths to be included in the 
evaluation of ^^(t). For a spin-isospin independent interaction the wave function is a scalar, 
and one can perform a fixed-node calculation in which configurations where the variational 
wave function is zero are discarded. This defines a surface within which the evolution 
proceeds, and eliminates the sign problem at the cost of introducing an approximation into 
the calculation. The fixed-node method is exact when the variational wave function has the 
true nodal surface, and provides an upper bound to the true ground state energy. This upper 
bound is often quite accurate because we are solving for the "optimum" solution subject 
only to the boundary condition that the wave function is zero on a predefined surface close 
to the correct one. 

The nuclear case is more complex, both because the trial wave function is a set of complex 
amplitudes and because we cannot evaluate the full wave function for a given set of coordi- 
nates. We can only evaluate it for a specific order of pair correlation operators in Eg 12. 4|, 
as a complete set would require [y4(y4 — l)/2]! terms. These pair orders are sampled in both 
the VMC and GFMC calculations. Fortunately the fluctuations in samples of pair orders 
arise from the commutators of correlation operators. These involve clusters of three or more 
nucleons, and they have a small effect on the variance. For the nuclear case, we construct an 
alternative constraint based on the overlap of each configuration \E'(r, Rj) with the sampled 
variational wave function. Configurations with negative overlap are discarded along with 
those with correspondingly small positive overlap, ensuring that the average overlap of the 
discarded configurations with the trial wave function is zero. This yields a stable simulation, 
and the calculation can proceed out to quite large imaginary time, much larger than the 
inverse gap in the system. This approximation is not guaranteed to produce an upper bound 
to the ground-state energy, though it has proven to be quite accurate for few-body nuclei. 
The results obtained using this method are labeled with CP. 

The CP approximation is tested by removing the constraint. The configurations gener- 
ated by the CP calculation are evolved further in the imaginary time r without constraint. 
The fermion sign problem makes this calculation more difficult for increasing density and 
for increasing r. In principle we can evaluate the energy for a much larger r at low density. 
In practice the low-density calculations appear to be well converged at fairly small imag- 
inary time. Of course the total unconstrained imaginary time propagation is quite small 
here, typically 0.005 MeV~^, and hence only fairly high-energy excitations are removed by 
this procedure. The results of these unconstrained (UC) GFMC calculations are the most 
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accurate of the presented results. 

In this work we use r to denote the time after the CP propagation. The CP GFMC 
propagation starts at a large negative r and ends at r = 0. The propagation time of CP 
GFMC is large enough to ensure convergence; however that of UC GFMC is limited by the 
growth in statistical errors due to the fermion sign problem. 

Results for the VMC and GFMC calculations at different densities are presented graph- 
ically in Figure [T] The upper square point at r = 0, at each density, is the VMC results 
for the v6 Hamiltonian, while the lower square point, also at r = 0, shows the CP GFMC 
results. The CP GFMC energies are typically 5 - 10% lower than the VMC for the v6 inter- 
action. The f 6 VMC and the PB variational chain summation (VCS) calculations discussed 
in section IV use the same wave function. 

The UC GFMC results are plotted as a function of unconstrained propagation time r 
after the end of CP propagation. Circles and squares show results for the v8' and v6. At all 
densities, the v6 calculations appear to be fairly stable and little change is observed between 
the CP results and the unconstrained results for larger imaginary time. Table H] lists the total 
and potential energy per neutron for various densities. The GFMC potential energies are 
approximately 15% lower than the VMC results indicating that true ground state has more 
correlations than the present variational wave function. The difference in the VMC and the 
GFMC potential energies is more than twice that in the total energies as expected near the 
minimum. The results of VCS calculations are also listed in Table U they are discussed in 
section IV. 

The results for v8' interaction are given in Table [H] and Figure ^ The VMC rows in this 
table give results with the variational wave function for the v6 potential without any spin- 
orbit correlation. With this wave function the expectation value of the spin-orbit interaction, 
(vl-s), is small and positive. It in nonzero due to the tensor correlations. In contrast the 
variational wave function used in the VCS calculations contain spin-orbit correlations which 
give significant negative (vl-s)- The L ■ S correlations absent in the v6 variational wave 
function are partly generated via the CP propagation as can be seen from the GFMC- 
CP (vl-s) values. However, the constraint imposed by the v6 wave function hinders their 
growth. After the constraint is removed, the spin-orbit correlations increase substantially, 
and we obtain significantly more attaction from the vl-s- This will be more evident in the 
comparison of pair distribution functions in section VI. The UC GFMC energy decreases 
with T (Fig. and at p > po the growth in statistical error limits the UC calculation. 

We have also performed calculations with different input correlation functions in the trial 
wave function. Their results for p = po are illustrated by the two sets of square points in 
Figure 121 The points labeled GFMC(LR) are obtained with the trial wave function having 
pair correlation functions of range L/2, while those labeled GFMC(SR) have much shorter 
range input pair correlation functions. There appears to be very little dependence of the 
UC GFMC results upon the choice of the range of input two-body correlation functions; this 
has been checked for the pair distribution functions as well. 

We have also implemented the auxiliary-field Diffusion Monte Carlo (AFDMC) method 
of Schmidt and Fantoni jl^. Since this method scales much better with A than the GFMC 
method discussed here, it can be used to treat much larger systems. At present the trial wave 
function used in these calculations includes only spin-independent Jastrow factors times a 
Fermi-Gas determinant: 

^j=[l[nr,j)]^FG. (3.3) 

i<j 
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Each configuration now has 14 two-component vectors describing the relative amphtude and 
phases of the spin of each neutron. These spins rotate in the presence of fluctuating fields 
which, when summed, reproduce exactly the results of the two-nucleon interaction. As in 
the GFMC calculation, a constraint is imposed requiring a positive overlap between the 
configuration at any time r and the trial wave function. 

The two sets of UC AFDMC results shown in Figure |21 are obtained with two different 
estimators of the ground-state energy. The growth energy (black dots) is determined from 
the rate of increase/decrease in the population with imaginary time r, while the mixed 
estimate (red dots) is determined by the overlap of the configurations with the Hamiltonian 
acting on the trial wave function. These two estimates should be equal within statistical 
errors for small values of the time step At. 

The \E'j is a very simple trial function and hence does not provide an accurate constraint. 
The CP AFDMC energies at r = are higher than the CP GFMC because of this relatively 
poor constraint. As r increases beyond the CP propagation region, the energy drops and 
becomes compatible with the GFMC results. The statistical errrors are somewhat worse, 
though, as each configuration contains only a single set of 14 spin vectors rather than the 2^^ 
amplitudes in the GFMC. The correlations between these amplitudes reduce the fermion sign 
problem, but at the cost of an exponentially increasing computational time. The AFDMC 
method has been used to study much larger systems with this simple constraint, and also to 
study the spin susceptibility of neutron matter. It could be used to determine the difference 
between the infinite-particle limit and the results for 14 neutrons. Here, though, we use 
VCS methods to calculate this difference. In addition, the QMC results provide a test of 
the VCS calculations often used in studies with more realistic Hamiltonians that include 
three-nucleon interactions and relativistic corrections. 



IV. VARIATIONAL CHAIN SUMMATION CALCULATIONS 

In VCS calculations of UG the expectation value of H — Tpc'- 

is expanded in powers of the short range functions {Fij — 1) |23|]. The ^fg is an eigenstate 
of the kinetic energy operator T = — V^/2m with the eigenvalue TpG = 0.3 k'^p/m, hence 
the terms with T operating on ^pg are not included in the expansion. The ra-body cluster 
contribution contains all the terms of this expansion having n neutrons. 

The leading two-body cluster contribution to the energy of UG of neutrons is given by: 

E{2h) = P-j d% C[F(r,,) (-i-V^ + v{n,)^ F(r,,)] 

- P^J d% C[e,, F(r,,) VF(r,,)] ■ i{r,j)Vi{n,) . (4.2) 

Here C[...] denotes the spin independent part, called the C-part 23] of the operators inside 
the square parenthesis, Cij is the spin exchange operator: 

= -1 (1 + <Ti ■ (Tj) , (4.3) 
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and i{r) is the spatial density matrix: 



^(r) 




,ikir 



(4.4) 



normalized such that i{r = 0) = 1. It is 



given by the Slater function: 



i{r) = 3[sin{x) 



X cos{x)]/x^ ; 



X = kpr . 



(4.5) 



for the $FG- 

It is relatively simple to calculate the above 2-body cluster contribution without approx- 
imations. All the terms in the 3-body cluster energy except those containing spin-orbit 
correlations can now be calculated exactly P|. However, all the > 4-body cluster contribu- 
tions as well as the 3-body contributions from L ■ S correlations are estimated approximately 
using the chain summation methods. 

Results of VCS calculations of the UG are given in Table ITTTl These are at optimum values 
of d and df, which generally exceed the L/2 of 14 neutron PB. In this case the Ey obtained 
with /3p = 1 is within ~ 2 % of the Ey with optimum f3p. The contributions of clusters 
are calculated following jOj]. The 3- and > 4-b-static contributions do not include spin-orbit 
interaction and correlation terms; their contributions are listed in row > 3-b-L-S. The values 
of 3-b-static contributions calculated with the chain summation approximation are also given 
in Table UTTl for comparison. They are typically within 10 % of the exact values. The listed 
values of > 4-b-static contributions include the elementary 4-body circular exchange diagram 
discussed in Appendix A. It was omitted in previous calculations because it is generally 
small in symmetric nuclear matter. However, this contribution contains the factor s~^, where 
s = 2,4 is the spin-isospin degeneracy factor in neutron and symmetric nuclear matter. It 
is relatively larger in neutron matter, and its values are listed in Table ITTTl 

Table ITTTl clearlv shows that the cluster expansion of the E{p) has slow convergence. At 
low densities this is primarily due to the large nn scattering length. With the v8' interaction, 
the the total contribution of clusters with > 4 is ~ 30 (10) % of the total energy at p = 



Table ITVl gives the results for UG variational energy for d = dt = L/2 and optimum 
values of a and f]p. The cluster expansion has better convergence for these shorter range 
correlations, and the Ev{d = dt = L/2) is above the optimum Ey by only ~ 0.1 MeV for 
P < PO) while at p = 1.5po it is higher by 0.3 MeV. At small p the j3t is significantly larger 
than 1. 

A. Variational Calculations with ^pb 

These calculations use the truncated interaction, v{rij)9{L/2—rij), and correlation ranges 
d = dt = L/2. The ^pb is an eigenstate of the kinetic energy with the eigenvalue (per 
neutron) : 



We note that the kinetic energy, Tpp per neutron of the 14 noninteracting neutrons in a 
periodic box is only 1.4 % larger than that of free Fermi gas. 



0.04 (0.16). 




(4.6) 
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As in the case of the UG we expand the expectation value of H — TpB as a sum over 
clusters. The leading two-body cluster contribution is given by Eq. ()4.2j) with the PB density 
matrix: 

^PB(r) = -(1 + 2 cos{xkB) + 2 cos{ykB) + 2 cos{zkB)) (4.7) 

in place of the UG density matrix given by the Slater function (Eq. 14. Sj) . The CpB depends 
upon the direction of r, as illustrated in Fig. El The Cpb is largest for r parallel to the box 
side and smallest along the diagonal. 

The C-parts of the operators in E{2b) (Eq. 14. 2j) are spherically symmetric functions of 
Tij, therefore the E{2h) depends only on the angle averaged value of 

ilB{.r) = ^ / sine de d(j) fpsir) , (4.8) 

4:71 J 

at r < L/2. The (-"ppi^) is fairly close to i^{r) in UG as shown in Fig. El Therefore 
the 2-body cluster contributions obtained with the UG and PB density matrices are not 
very different. In the variational calculations with $pb we approximate the contributions 
of all n > 3-body clusters by their values in the UG. The 3-b-static cluster is calculated 
exactly, and the rest with chain summation approximations. We note that Fantoni and 
Schmidt [24!\ have developed chain summation methods for calculations with ^pb sans tensor 
correlations. They retain the ipB in all the many-body cluster contributions calculated with 
chain summation methods. 

The results of calculations with the v6 and v8' interactions are given in Tables IVl and IVll 
The values of Tpc — TpB and 

AE{2b) = Euci^b) - EpB{2h) (4.9) 

are also listed in Table IVTl The smallness of these differences makes the 14- neutron periodic 
box a good approximation for studying uniform gas E{p). The variational parameters a and 
j3p have essentially the same values in PB calculations as in UG with d = dt = L/2. The main 
difference between the UG and PB energies comes from the contribution of the long range 
interaction, v{rij)6{rij — L/2) omitted in the PB. Its contribution denoted by {vijij > L/2)) 
is estimated from UG calculations and listed in Table IVII It becomes comparable to the 
total Eug{p) at p ~ 1.5 po- 

B. Comparison with QMC Calculations 

The results of the approximate VCS calculations are compared with those of QMC cal- 
culations in Tables U and O for the v6 and v8' interactions. Since the VMC and VCS 
calculations use the same wave function for the v6 Hamiltonian, they should ideally give the 
same results. The difference between them is due only to the approximations in the VCS 
calculations. Relative to the VMC results, the VCS total energy is higher by 8 % at po/4 
and lower by 2 % at 1.5po. The difference between the VCS and VMC potential energy, (w6) 
is similar. The clusters with > 3-neutrons give a relatively larger contribution at smaller 
densities (Table [VJ due to the large nn scattering length. It is thus not surprising that VCS 
has larger errors in that region. The GFMC-UC energies are below those of the VCS by 
12 to 4 % in this density range. We note that the differences between the GFMC-UC and 
VCS or VMC potential energies are much larger, of order 20 %. This is because the present 
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variational wave functions underestimate the correlations in matter, as is further elaborated 
in the next section. 

In the case of the v8' Hamiltonian, the VMC calculations do not include the spin-orbit 
correlations while the VCS do, therefore the VMC and VCS results are not comparable 
in this case. However, we can compare the GFMC-UC and VCS results. As for v6, the 
GFMC-UC energies are lower by 10 % at lower densities; however, at 1.5po the VCS energy 
is lower by 2±0.5 MeV. Most of this difference seems to come from the spin-orbit interaction. 
The (fL-s) = — 8 ± 1 and —12 MeV in these two calculations at l.Spo- When the density 
dependence of spin-orbit correlations is neglected, the leading two-neutron cluster gives 
a contribution proportional to p^/^ to (fL s)- This comes about because the spin-orbit 
contributions are proportional to p^^^ via the fc^ momentum dependence of wl sL ■ S/l-sL ■ S, 
and the summation over particles gives an additional factor of p for 2-body clusters. The 
VCS results for (vl s) approximately follow this density dependence as shown in Table ITTl 
Up to po the GFMC-UC (fL s) also has a similar density dependence. However, at 1.5po 
the GFMC-UC is smaller in magnitude. It could be that higher order cluster terms become 
more important at this density and that VCS overestimates the wl-s contribution, or that, 
if GFMC-UC is propagated further in the imaginary time r after the CP, the (fLs) will 
decrease and the GFMC-UC energy will go down. In the present calculation we can not 
test this possibility because of the increase in the statistical errors due to the Fermion sign 
problem. 



V. PAIR DISTRIBUTION FUNCTIONS 

The pair distribution functions obtained from QMC calculations with the v8' interaction 
are shown in Figures |3] to [7| In each figure the circles show the results of VMC calculations 
with a wave function containing /6 correlations (without L-S correlations). Squares and 
triangles represent the result of the constrained path (GFMC-CP), and the unconstrained 
(GFMC-UC). The pair distribution functions gp{r) are given by the expectation value: 

9pir) = NY^{<il\5in, - r)Or,.|vP), (5.1) 

with a normalization such that gi{r) = gdr) goes to one at large distances. The (72-4 are 
denoted by ga,gt and g^s for clarity. The gdr) gives the probability to find a neutron at 
a distance r from another neutron since Of^"^ = 1. In contrast Of~^ = CTi ■ (Tj, thus g^ is 
proportional to the expectation value of 6{rij — r)(Ti ■ CTj. Using the projection operators: 



^^5=0 = \{l-(T,-(Tj) , (5.2) 



PS=1 = \{?> + <Ti-(T,) , (5.3) 

the pair distribution functions in spin 5 = and 1 pairs are found to be: 

gs=Q{r) = ^{gc{r) - ga{r)) , (5.4) 

gs=i{r) = ^{^gc{r) + ga{r)) . (5.5) 
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Since gs=i{f^ ~^ 0) — 0, the gcr{r) = —3gc{r) at small r. In noninteracting Fermi gases: 

g^^ir) = 1 - ^fir) , (5.6) 
g^'^ir) = —fir) , (5.7) 

g^ir) = g[f = . (5.8) 

The VMC calculations do not have spin-orbit correlations, and give gis ~ 0. 

At p = 0.04 fm~^, there is a very strong pairing into spin states as indicated by the 
large negative g^- It can be seen more clearly in Figure |H1 which compares the gs=o,i{f') in 
neutron matter and Fermi gases. The large peak of the gs=o is due to the large negative nn 
scattering length; it should be relatively model independent and grow at smaller densities. 
This pairing is present in the variational calculations, though underestimated by ~ 25 %. 
The tensor and L-S correlations are quite modest at these low densities. There is little change 
between the constrained GFMC-CP and unconstrained GFMC-UC, indicating reasonably 
good convergence within this class of wave functions. The tensor correlations are long- 
ranged, extending nearly to the L/2 limit imposed by the periodic boundary conditions. 
This same behavior is seen even when starting the GFMC with trial wave functions having 
much shorter range correlations. 

The correlations at p = 0.08 fm~'^ are fairly similar, though the spin correlation is not 
as large, and the tensor and L ■ S correlations are becoming more significant. The gLs is 
essentially zero in the variational calculation, and underestimated in CP GFMC. 

At the largest densities considered, p = 0.16 and 0.24 fm~^, the differences between the 
variational, GFMC-CP, and the GFMC-UC results are quite large. Both the tensor and L-S 
correlations are quite important and significantly underestimated in the VMC and GFMC- 
CP calculations. We see a transition from low densities, where the S-wave interaction and 
spin zero pairing is dominant, to these higher densities, where the P-wave interactions are 
crucial. It could be associated with the ^P2-^F2 pairing [i^ expected at higher densities. 
In VCS calculations with three nucleon interaction (see Fig. 9 of Ref . |26l| ) such a behavior is 
associated with the onset of pion condensation. 



VI. DENSITY DEPENDENCE OF NEUTRON MATTER ENERGY 

The total energy of neutron matter interacting with the v8' potential is reported in Table 
I VI II It is obtained by adding the box corrections listed in Table IVII to the GFMC-UC 
energies listed in Table ITTl The ratio of neutron matter E{p) to the noninteracting neutron 
Fermi gas energy is also listed in Table Ivnl This ratio approaches ~ 0.5 at low densities. 

The properties of low density neutron matter are dominated by the large negative nn 
scattering length. When \akp\ « 1 we have the well known low density expansion [27l |: 



E{p) = Efg{p) 



10 4 
1 + —akp + 77-^(11 - 2ln2){akFf + ... 



(6.1) 



Such an expansion is not useful for neutron matter because even at densities as low as 1% 
of po, \akp\ > 6. The limit akp — oo is perhaps more applicable to neutron gas than the 
low density expansion, as suggested by Bertsch [15|. In this limit it is known that: 

E{p) = EpG{p)i ■ (6.2) 
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The estimates of ^ range from 0.326 H to 0.568 to 0.59 Recent quantum 

Monte Carlo calculations give ^ ~ 0.54 for the normal phase and ^ = 0.44 ± 0.02 for 
the superfluid phase. Most many-body calculations, both Brueckner and variational, give 
^ ~ 0.5 for normal neutron matter. As an example, we compare the energies of neutron 
matter calculated with the CSM in 1981 with the Epc and results of present calculation 
in Figure ini 

Eq. 16.21 implies that at low densities the interaction energy of neutron matter becomes 
proportional to kp as is the FG kinetic energy. This interaction energy is proportional to 
density (= kp/Sn^) times the volume integral of the effective G-interaction, related to the 
bare f-interaction by the well-known Brueckner equation G(j) = vip. Here and ip are 
the unperturbed and perturbed two-nucleon wave functions. At small relative momenta of 
interest at low-densities, = 1, and in vacuum = 1 — a/r beyond the the range of 
V. When —a/R^ » 1, as is the case for neutrons, we can neglect the 1 in comparison and 
approximate the ip by —a/R^. The effective interaction in vacuum is essentially enhanced 
by a factor —a/Ry by the large scattering length. In matter the effective scattering length 
is limited by the interparticle spacing of order l/kp. Thus, when —akp » 1, the G is 
enhanced by a factor proportional to l/kpRv, its integral becomes proportional to l/kp, 
and the interaction energy proportional to kj^. At higher densities we see a deviation from 
Eq. 16.21 in Fig. IHl It starts when kpRv becomes of order 1 and the first (unit) term of ip can 
not be neglected. When i?^ — >^ 0, as in the challenge problem proposed by Bertsch, Eq. 16.21 
is valid at all densities when a = — oo. 

Most of the nonrelativistic Skyrme as well as the relativistic mean field energy density 
functions commonly used to study nuclei and neutron star matter assume that ^ = 1. None 
of these therefore can reproduce the equation of state, Eq{p) of pure neutron matter [l| 
obtained from realistic interactions even at low densities. The effective interaction used in 
mean field models must diverge as p — > due to large wi scattering length. Energy density 
functionals containing such low-density divergences [31] are probably necessary to study 
nuclei near neutron drip line or in the inner crust of neutron stars. 



VII. SUPERFLUIDITY AND ACCURACY OF THE PRESENT CALCULATIONS 

Quantum Monte Carlo calculations of 14 neutrons described above start from a correlated 
Slater trial wave function. This trial function is appropriate for the normal phase of Fermi 
liquids and is used here as the starting point for the CP GFMC and to obtain the constraints 
used to limit the Fermion sign problem in that calculations. Its results are expected to give 
the equation of state of the normal phase. The long-wavelength properties of the correlated 
Slater wave function do not include the expected superfluid properties of neutron matter. 
Generally, the energy of the superfluid phase is not significantly different from that of the 
normal one, since in most systems the pairing is confined near the Fermi surface, and involves 
only relatively few particles. Here, however, the magnitude of the nn scattering length is 
very large compared to the interparticle spacing, the pairing is exceptionally strong, and 
affects all the particles. 

In principle it may be possible for the unconstrained GFMC calculations to relax to 
the lower energy superfluid phase provided there is sufficient overlap between the fourteen 
particle wavefunctions of these phases. However, the unconstrained calculations are limited 
to fairly short paths, due to small unconstrained propogation time, and hence are unlikely 
to relax to the superfluid phase. In calculations with simple spin-independent, short-range 
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interactions with a = — oo, Schmidt et al. [17| find, upon including superfluid pairing into 
the trial wave function, important effects, including an ^ 20% reduction in the energy per 
particle mentioned in the previous section. 

One of the important features of the Monte Carlo approach pursued here is that it can 
be extended to include BCS pairing into the trial wave function. Thus it will be possible 
to study properties of superfluid neutron matter with realistic models of the nn interaction. 
We are presently pursuing such studies. The main remaining uncertainty in the equation of 
state of low density neutron matter is probably due to the neglected difference in the Eo{p) 
of normal and superfluid phases. 
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APPENDIX A: ELEMENTARY FOUR-BODY CIRCULAR EXCHANGE DIA- 
GRAMS 

The 4-body diagrams in the expansion of Ey — Tpc (Eq. 14. 1|) . not included in the Fermi 
hypernetted chain summation, are called elementary [23] . In Fermi liquids they were first 
calculated by Zabolitzky js^]. The circular exchange diagram shown in Fig. is the 
only one of first order in the expansion in powers of (F,^j — 1), among these. It therefore 
has special importance as emphasized by Krotscheck j33|, and was included in the study of 
nuclear matter structure functions [s^. The thick interaction line in this diagram represents: 




The term with 634623612 takes into account the circular exchange in the other direction. The 
spin-orbit interactions and correlations, and the spin correlations between neutron pairs 12, 
23, 34 and 41 are neglected in this approximation expected to have an accuracy of order 20 





(Al) 



We approximate the contribution of this diagram with: 




(A2) 



% 
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TABLE I: Quantum Monte Carlo and VCS results for 14 neutrons in PB with the v6 Hamiltonian 
(MeV per neutron). Statistical errors are indicated in parenthesis. 



Method 


p (fm ^) 


0.04 


0.08 


0.16 


0.24 


VMC 


{H) 


7.04(01) 


11.32(01) 


21.39(01) 


34.30(01) 


GFMC-CP 




6.72(01) 


10.64(01) 


19.80(02) 


31.90(02) 


GFMC-UC 




6.75(01) 


10.64(03) 


19.91(11) 


32.15(08) 


VCS 




7.6 


11.9 


21.2 


33.6 


VMC 




-9.92(03) 


-16.17(05) 


-22.79(08) 


-26.10(11) 


GFMC-CP 




-11.75(08) 


-18.64(09) 


-28.01(08) 


-32.94(25) 


GFMC-UC 




-11.36(13) 


-18.17(36) 


-26.62(71) 


-32.71(71) 


VCS 




-9.2 


-15.4 


-22.7 


-26.6 



TABLE II: Quantum Monte Carlo and VCS results for 14 neutrons in PB with the v8' Hamiltonian 
(MeV per neutron). 



Method 


p (fm ^) 


0.04 


0.08 


0.16 


0.24 


VMC 


(H) 


7.16(01) 


11.678(07) 


21.82(12) 


35.02(01) 


GFMC-CP 




6.43(01) 


10.02(02) 


18.54(04) 


30.04(04) 


GFMC-UC 




6.32(03) 


9.501(06) 


17.00(27) 


28.35(50) 


VCS 




7.0 


10.3 


17.4 


26.3 


VMC 


{v6) 


-9.74(03) 


-15.72(6) 


-21.64(09) 


-24.37(11) 


GFMC-CP 




-11.85(09) 


-18.34(11) 


-27.72(15) 


-32.34(23) 


GFMC-UC 




-11.44(19) 


-17.83(30) 


-25.62(87) 


-30.52(1.35) 


VCS 




-9.3 


-15.1 


-21.6 


-24.9 


VMC 


{VL-S) 


0.07(01) 


0.26(01) 


0.13(01) 


0.21(01) 


GFMC-CP 




-0.23(11) 


-1.37(03) 


-2.69(03) 


-4.08(08) 


GFMC-UC 




-0.85(04) 


-2.59(15) 


-6.24(50) 


-7.98(98) 


VCS 




-0.88 


-2.3 


-6.9 


-12.1 


-6.9(p/po)'/' 




-0.68 


-2.2 


-6.9 


-13.6 
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TABLE III: Results of VCS calculations of the UG E{p) (MeV per neutron) with optimum d 
and dt, in MeV. Contributions with an * are estimated using chain summation approximation; 
those without are calculated exactly. The bottom two rows give the approximate 3-b-static*, for 
comparison with the exact 3-b-static, and 4-b-elementary* circular exchange contribution included 
in > 4-b-static* 



p (fm '*) 


0.04 


0.08 


0.16 


0.24 


a (imj 


6.00 


o.ol 


z.z9 


O OA 

2.20 


at (rmj 


0.1 ( 


0.00 


o.zz 


4. by 


a 


U.oy 


v.oi 


n on 


0.(2 


Pa 


U.o 


U.o 


i.O 


i.O 


A 


0.9 


1.0 


0.9 


0.9 




0.9 


0.8 


1.0 


1.1 




13.9 


22.1 


35.1 


46.0 


2-b-Total 


-10.1 


-16.9 


-25.5 


-33.9 


3-b-static 


4.9 


6.9 


3.7 


3.5 


> 4-b-static* 


-2.3 


-3.1 


- 1.5 


- 1.1 


> 3-b-L • S* 


0.1 


0.2 


0.2 


- 0.1 


Total E{p) 


6.6 


9.2 


12.0 


14.5 


3-b-static* 


4.6 


6.6 


3.5 


3.2 


4-b-elementary* 


- 0.5 


-0.5 


0.5 


0.7 



TABLE IV: Results of VCS calculations of the UG E{p) (MeV per neutron) with d = dt = L/2, in 
MeV. Contributions with an * are estimated using chain summation approximation; those without 
are calculated exactly. 



p (fm-3) 


0.04 


0.08 


0.16 


0.24 


L/2 (fm) 


3.52 


2.80 


2.22 


1.94 


a 


0.90 


0.85 


0.80 


0.80 


Pa 


0.80 


0.9 


0.9 


1.0 


A 


1.50 


2.0 


1.0 


1.0 




0.85 


0.9 


1.1 


1.1 


TpG 


13.9 


22.1 


35.1 


46.0 


2-b-Total 


- 9.7 


-14.9 


-23.0 


-29.9 


3-b-static 


3.9 


2.7 


-0.1 


-0.8 


> 4-b-static* 


-1.5 


-0.7 


- 0.0 


- 0.0 


> 3-b-L • S* 


0.0 


0.0 


0.2 


- 0.4 


Total E{p) 


6.7 


9.2 


12.1 


14.8 
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TABLE V: Results of PB variational calculations with v6 interaction truncated at r = L/2. The 
variational parameters are listed in Table ITVl and the energies are MeV per neutron. 



p (fm ^) 


0.04 


0.08 


0.16 


0.24 


TpB 


14.1 


22.4 


35.6 


46.6 


2-b-Total 


-9.0 


-12.6 


-14.2 


-12.3 


>3-b-Total 


2.5 


2.0 


-0.1 


-0.7 


Total E{p) 


7.6 


11.9 


21.2 


33.6 



TABLE VI: Results of PB variational calculations with u8' interaction truncated at r = L/2. The 
variational parameters are listed in Table HVl and the last three rows give the differences between 
PB and UG contributions. All energies are in MeV per neutron. 



p (fm ^) 


0.04 


0.08 


0.16 


0.24 


TpB 


14.1 


22.4 


35.6 


46.6 


2-b-Total 


-9.5 


-14.0 


-18.3 


-19.1 


>3-b-Total 


2.4 


1.9 


0.1 


-1.2 


Total E{p) 


7.0 


10.3 


17.4 


26.3 


TpG — TpB 


-0.2 


-0.3 


-0.5 


-0.7 


AE{2b) 


-0.0 


-0.1 


-0.1 


-0.1 


{v{n, > L/2)) 


-0.1 


-0.8 


-4.5 


-10.7 



TABLE VII: Neutron matter energy with the w8' interaction in MeV per neutron. 



p (fm '^) 


0.04 


0.08 


0.16 


0.24 


GFMC-UC 


6.3 


9.5 


17.0 


28.4 


Box Correction 


-0.3 


-1.1 


-5.1 


-11.5 


Total E{p) 


6.0 


8.4 


12.1 


16.9 


E{p)/Efg{p) 


0.43 


0.38 


0.34 


0.37 
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FIG. 1: Energy vs. imaginary time r after CP propagation, at various densities. VMC results 
(upper squares) for v6 and GFMC-CP results (lower squares and dots) for vG and v8' are shown 
at r = 0, and unconstrained GFMC results are shown for various r > 0. 
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FIG. 2: Energy vs. imaginary time r for different calculations using v8' interaction at p = 0.16 
fm^"^. Two different estimates (growth and mixed) are shown for the AFDMC calculation along 
with the results of two different GFMC calculations, using short-range (SR) and long-range (LR) 
correlations in the trial wave function. 
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r (fm) 

FIG. 3: The Slater function lpB{r) for 14 neutrons in PB. The top and bottom dash-dot lines 
show the •^ps(r) with r parallel to the 0,0,1 box side and to the 1,1,1 diagonal respectively. The 

middle full line shows the average v/^p£(r), and the dashed line shows the i{r) in the UG. 
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FIG. 4: Pair distribution functions for the v8' interaction at rho = 0.04 fm^'^. The data sets from 
top to bottom correspond to central, L-S, tensor and <Tj • crj pair distribution functions. Each set 
contains circles, squares and triangles showing the VMC, GFMC-CP and GFMC-UC results. The 
L-S and the tensor distribution functions are scaled up by a factor of five. 



E 

Q. 



1.0 
0.5 
0.0 
-0.5 
-1.0 



i 







g g s ^ ^ ^ 

^ H H 





^ 







■"■5^ VMC 

□ GFMC-CP 
-2.0 GFMC-UC 

"^■%.0 ' o!5~ 



® © © ^ ^ S 

A ^ 7X ^ ^ 

11 S ^ lU 

fi B 



_L 



_L 



1.0 1.5 

r(fm) 



G 

mi 






9i 


@ 


9a 





g, (x5) 


© 


gLs (x5) 



AV8' 

p = 0.08 fm' 



2.0 



2.5 



FIG. 5: Pair distribution functions for the v8' interaction at rho = 0.08 fm ^, as in Figure [l] 
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FIG. 6: Pair distribution functions for the v8' interaction at rho = 0.16 fm ^, as in Figure^ 
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FIG. 7: Pair distribution functions for the v8' interaction at rho = 0.24 fm as in Figure H) 
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FIG. 8: Pair distribution functions for spin and spin 1 pairs at p = 0.04 fm~^; results of 
unconstrained GFMC calculations are compared to distributions in noninteracting FG shown by 
solid lines. 
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FIG. 9: The energy of neutron matter at low densities. The stars give the results of Ref.j30|] and 
the circles of the present calculation. The dashed line shows 0.5 Epc- 
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FIG. 10: The lowest order elementary 4-body circular exchange diagram. The dashed line shows 
F24 — 1, the thick solid line denotes the interaction link, and the thin lines with direction show the 
exchange pattern. 
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